home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / TPL60N14.ARJ / MACHARIT.PAS < prev    next >
Pascal/Delphi Source File  |  1992-05-01  |  14KB  |  368 lines
  1. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  2.  
  3. UNIT MachArit;   { converted from Fortran original 05-01-92 Norbert Juffa }
  4.  
  5. INTERFACE
  6.  
  7. PROCEDURE MACHAR (VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  8.                       MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
  9.  
  10. PROCEDURE PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  11.                       MAXEXP: LONGINT; EPS,EPSNEG,XMIN,XMAX: REAL);
  12.  
  13.  
  14. IMPLEMENTATION
  15.  
  16. PROCEDURE MACHAR(VAR IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  17.                      MAXEXP: LONGINT; VAR EPS,EPSNEG,XMIN,XMAX: REAL);
  18.  
  19.  
  20. {-----------------------------------------------------------------------
  21. C  This Fortran 77 subroutine is intended to determine the parameters
  22. C   of the floating-point arithmetic system specified below.  The
  23. C   determination of the first three uses an extension of an algorithm
  24. C   due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
  25. C   but not all, of the improvements suggested by M. Gentleman and S.
  26. C   Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
  27. C   program was published in the book Software Manual for the
  28. C   Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
  29. C   Englewood Cliffs, NJ, 1980.  The present version is documented in
  30. C   W. J. Cody, "MACHAR: A subroutine to dynamically determine machine
  31. C   parameters," TOMS 14, December, 1988.
  32. C
  33. C  The program as given here must be modified before compiling.  If
  34. C   a single (double) precision version is desired, change all
  35. C   occurrences of CS (CD) in columns 1 and 2 to blanks.
  36. C
  37. C  Parameter values reported are as follows:
  38. C
  39. C       IBETA   - the radix for the floating-point representation
  40. C       IT      - the number of base IBETA digits in the floating-point
  41. C                 significand
  42. C       IRND    - 0 if floating-point addition chops
  43. C                 1 if floating-point addition rounds, but not in the
  44. C                   IEEE style
  45. C                 2 if floating-point addition rounds in the IEEE style
  46. C                 3 if floating-point addition chops, and there is
  47. C                   partial underflow
  48. C                 4 if floating-point addition rounds, but not in the
  49. C                   IEEE style, and there is partial underflow
  50. C                 5 if floating-point addition rounds in the IEEE style,
  51. C                   and there is partial underflow
  52. C       NGRD    - the number of guard digits for multiplication with
  53. C                 truncating arithmetic.  It is
  54. C                 0 if floating-point arithmetic rounds, or if it
  55. C                   truncates and only  IT  base  IBETA digits
  56. C                   participate in the post-normalization shift of the
  57. C                   floating-point significand in multiplication;
  58. C                 1 if floating-point arithmetic truncates and more
  59. C                   than  IT  base  IBETA  digits participate in the
  60. C                   post-normalization shift of the floating-point
  61. C                   significand in multiplication.
  62. C       MACHEP  - the largest negative integer such that
  63. C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
  64. C                 MACHEP is bounded below by  -(IT+3)
  65. C       NEGEPS  - the largest negative integer such that
  66. C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
  67. C                 NEGEPS is bounded below by  -(IT+3)
  68. C       IEXP    - the number of bits (decimal places if IBETA = 10)
  69. C                 reserved for the representation of the exponent
  70. C                 (including the bias or sign) of a floating-point
  71. C                 number
  72. C       MINEXP  - the largest in magnitude negative integer such that
  73. C                 FLOAT(IBETA)**MINEXP is positive and normalized
  74. C       MAXEXP  - the smallest positive power of  BETA  that overflows
  75. C       EPS     - the smallest positive floating-point number such
  76. C                 that  1.0+EPS .NE. 1.0. In particular, if either
  77. C                 IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  78. C                 Otherwise,  EPS = (FLOAT(IBETA)**MACHEP)/2
  79. C       EPSNEG  - A small positive floating-point number such that
  80. C                 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
  81. C                 or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  82. C                 Otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
  83. C                 NEGEPS is bounded below by -(IT+3), EPSNEG may not
  84. C                 be the smallest number that can alter 1.0 by
  85. C                 subtraction.
  86. C       XMIN    - the smallest non-vanishing normalized floating-point
  87. C                 power of the radix, i.e.,  XMIN = FLOAT(IBETA)**MINEXP
  88. C       XMAX    - the largest finite floating-point number.  In
  89. C                 particular  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  90. C                 Note - on some machines  XMAX  will be only the
  91. C                 second, or perhaps third, largest number, being
  92. C                 too small by 1 or 2 units in the last digit of
  93. C                 the significand.
  94. C
  95. C     Latest revision - December 4, 1987
  96. C
  97. C     Author - W. J. Cody
  98. C              Argonne National Laboratory
  99. C
  100. C-----------------------------------------------------------------------}
  101.  
  102. VAR   I, L, ITEMP, IZ,
  103.       J, K, MX, NXRES:     LONGINT;
  104.       A, B, BETA, BETAIN,
  105.       BETAH, ONE, T, TEMP,
  106.       TEMPA, TEMP1, TWO,
  107.       Y, Z, ZERO:          REAL;
  108.       CONV:                ARRAY [0..10] OF REAL;
  109.  
  110. LABEL 1, 2, 10, 21, 22, 30, 32, 40,
  111.       41, 42, 43, 44, 45, 46, 50, 52;
  112.  
  113. BEGIN
  114.  
  115. {-----------------------------------------------------------------------}
  116.  
  117.    FOR L := 1 TO 10 DO BEGIN
  118.       CONV [L] := L;
  119.    END;
  120.  
  121.    ONE  := CONV [1];
  122.    TWO  := ONE + ONE;
  123.    ZERO := ONE - ONE;
  124.  
  125. {-----------------------------------------------------------------------}
  126. {  Determine IBETA, BETA ala Malcolm.                                   }
  127. {-----------------------------------------------------------------------}
  128.  
  129.    A := ONE;
  130. 1: A := A + A;
  131.    TEMP  := A+ONE;
  132.    TEMP1 := TEMP-A;
  133.    IF TEMP1-ONE = ZERO THEN
  134.       GOTO 1;
  135.    B := ONE;
  136. 2: B := B + B;
  137.    TEMP := A+B;
  138.    ITEMP := TRUNC (TEMP-A);
  139.    IF ITEMP = 0 THEN
  140.       GOTO 2;
  141.    IBETA := ITEMP;
  142.    BETA  := CONV [IBETA];
  143.  
  144. {-----------------------------------------------------------------------}
  145. {  Determine IT, IRND.                                                  }
  146. {-----------------------------------------------------------------------}
  147.  
  148.    IT := 0;
  149.    B  := ONE;
  150. 10:IT := IT + 1;
  151.    B := B * BETA;
  152.    TEMP := B+ONE;
  153.    TEMP1 := TEMP-B;
  154.    IF TEMP1-ONE = ZERO THEN
  155.       GOTO 10;
  156.    IRND := 0;
  157.    BETAH := BETA / TWO;
  158.    TEMP := A+BETAH;
  159.    IF TEMP-A <> ZERO THEN
  160.       IRND := 1;
  161.    TEMPA := A + BETA;
  162.    TEMP := TEMPA+BETAH;
  163.    IF (IRND = 0) AND (TEMP-TEMPA <> ZERO) THEN
  164.       IRND := 2;
  165.  
  166. {-----------------------------------------------------------------------}
  167. {  Determine NEGEP, EPSNEG.                                             }
  168. {-----------------------------------------------------------------------}
  169.  
  170.    NEGEP := IT + 3;
  171.    BETAIN := ONE / BETA;
  172.    A := ONE;
  173.    FOR I := 1 TO NEGEP DO BEGIN
  174.       A := A * BETAIN;
  175.    END;
  176.    B := A;
  177. 21:TEMP := ONE-A;
  178.    IF TEMP-ONE <> ZERO THEN
  179.       GOTO 22;
  180.    A := A * BETA;
  181.    NEGEP := NEGEP - 1;
  182.    GOTO 21;
  183. 22:NEGEP := -NEGEP;
  184.    EPSNEG := A;
  185.  
  186. {-----------------------------------------------------------------------}
  187. {  Determine MACHEP, EPS.                                               }
  188. {-----------------------------------------------------------------------}
  189.  
  190.    MACHEP := -IT - 3;
  191.    A := B;
  192. 30:TEMP := ONE+A;
  193.    IF TEMP-ONE <> ZERO THEN
  194.       GOTO 32;
  195.    A := A * BETA;
  196.    MACHEP := MACHEP + 1;
  197.    GOTO 30;
  198. 32:EPS := A;
  199.  
  200. {-----------------------------------------------------------------------}
  201. {  Determine NGRD.                                                      }
  202. {-----------------------------------------------------------------------}
  203.  
  204.    NGRD := 0;
  205.    TEMP := ONE+EPS;
  206.    IF (IRND = 0) AND (TEMP*ONE-ONE <> ZERO) THEN
  207.       NGRD := 1;
  208.  
  209. {-----------------------------------------------------------------------
  210. C  Determine IEXP, MINEXP, XMIN.
  211. C
  212. C  Loop to determine largest I and K = 2**I such that
  213. C         (1/BETA) ** (2**(I))
  214. C  does not underflow.
  215. C  Exit from loop is signaled by an underflow.
  216. C-----------------------------------------------------------------------}
  217.  
  218.    I := 0;
  219.    K := 1;
  220.    Z := BETAIN;
  221.    T := ONE + EPS;
  222.    NXRES := 0;
  223. 40:Y := Z;
  224.    Z := Y * Y;
  225.  
  226. {-----------------------------------------------------------------------}
  227. {  Check for underflow here.                                            }
  228. {-----------------------------------------------------------------------}
  229.  
  230.    A := Z * ONE;
  231.    TEMP := Z * T;
  232.    IF ((A+A = ZERO) OR (ABS(Z) >= Y)) THEN
  233.       GOTO 41;
  234.    TEMP1 := TEMP * BETAIN;
  235.    IF TEMP1*BETA = Z THEN
  236.       GOTO 41;
  237.    I := I + 1;
  238.    K := K + K;
  239.    GOTO 40;
  240. 41:IF IBETA = 10 THEN
  241.       GOTO 42;
  242.    IEXP := I + 1;
  243.    MX := K + K;
  244.    GOTO 45;
  245.  
  246. {-----------------------------------------------------------------------}
  247. {  This segment is for decimal machines only.                           }
  248. {-----------------------------------------------------------------------}
  249.  
  250. 42:IEXP := 2;
  251.    IZ := IBETA;
  252. 43:IF K < IZ THEN
  253.       GOTO 44;
  254.    IZ := IZ * IBETA;
  255.    IEXP := IEXP + 1;
  256.    GOTO 43;
  257. 44:MX := IZ + IZ - 1;
  258.  
  259. {-----------------------------------------------------------------------}
  260. {  Loop to determine MINEXP, XMIN.                                      }
  261. {  Exit from loop is signaled by an underflow.                          }
  262. {-----------------------------------------------------------------------}
  263.  
  264. 45:XMIN := Y;
  265.    Y := Y * BETAIN;
  266.  
  267. {-----------------------------------------------------------------------}
  268. {  Check for underflow here.                                            }
  269. {-----------------------------------------------------------------------}
  270.  
  271.    A    := Y * ONE;
  272.    TEMP := Y * T;
  273.    IF (A+A = ZERO) OR (ABS(Y) >= XMIN) THEN
  274.       GOTO 46;
  275.    K := K + 1;
  276.    TEMP1 := TEMP * BETAIN;
  277.    IF (TEMP1*BETA <> Y) OR (TEMP = Y) THEN
  278.       GOTO 45
  279.    ELSE BEGIN
  280.       NXRES := 3;
  281.       XMIN  := Y;
  282.       END;
  283. 46:MINEXP := -K;
  284.  
  285. {-----------------------------------------------------------------------}
  286. {  Determine MAXEXP, XMAX.                                              }
  287. {-----------------------------------------------------------------------}
  288.  
  289.    IF (MX > K+K-3) OR (IBETA = 10) THEN
  290.       GOTO 50;
  291.    MX := MX + MX;
  292.    IEXP := IEXP + 1;
  293. 50:MAXEXP := MX + MINEXP;
  294.  
  295. {-----------------------------------------------------------------}
  296. {  Adjust IRND to reflect partial underflow.                      }
  297. {-----------------------------------------------------------------}
  298.  
  299.    IRND := IRND + NXRES;
  300.  
  301. {-----------------------------------------------------------------}
  302. {  Adjust for IEEE-style machines.                                }
  303. {-----------------------------------------------------------------}
  304.  
  305.    IF IRND >= 2 THEN
  306.       MAXEXP := MAXEXP - 2;
  307.  
  308. {-----------------------------------------------------------------}
  309. {  Adjust for machines with implicit leading bit in binary        }
  310. {  significand, and machines with radix point at extreme          }
  311. {  right of significand.                                          }
  312. {-----------------------------------------------------------------}
  313.  
  314.    I := MAXEXP + MINEXP;
  315.    IF (IBETA = 2) AND (I = 0) THEN
  316.       MAXEXP := MAXEXP - 1;
  317.    IF I > 20 THEN
  318.       MAXEXP := MAXEXP - 1;
  319.    IF A <> Y THEN
  320.       MAXEXP := MAXEXP - 2;
  321.    XMAX := ONE - EPSNEG;
  322.    IF XMAX*ONE <> XMAX THEN
  323.       XMAX := ONE - BETA * EPSNEG;
  324.    XMAX := XMAX / (BETA * BETA * BETA * XMIN);
  325.    I := MAXEXP + MINEXP + 3;
  326.    IF I <= 0 THEN
  327.       GOTO 52;
  328.    FOR L := 1 TO I DO BEGIN
  329.       IF IBETA = 2 THEN
  330.          XMAX := XMAX + XMAX;
  331.       IF IBETA <> 2 THEN
  332.          XMAX := XMAX * BETA;
  333.    END;
  334. 52:
  335. END; { MachAr }
  336.  
  337.  
  338.  
  339. PROCEDURE PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  340.                       MAXEXP: LONGINT; EPS,EPSNEG,XMIN,XMAX: REAL);
  341. BEGIN
  342.    WRITELN ('MACHAR DETERMINED THE FOLLOWING PARAMETERS OF THE FLOATING-POINT ARITHMETIC');
  343.    WRITELN;
  344.    WRITELN ('RADIX FOR FLOATING-POINT REPRESENTATION:        ', IBETA);
  345.    WRITELN ('NUMBER OF BASE ', IBETA:2, ' DIGITS IN SIGNIFICAND:        ', IT);
  346.    IF IBETA <> 10 THEN
  347.       WRITELN ('NUMBER OF BITS USED TO REPRESENT THE EXPONENT:  ', IEXP)
  348.    ELSE
  349.       WRITELN ('NUMBER OF DECIMAL PLACES USED FOR THE EXPONENT:  ', IEXP);
  350.    WRITELN ('SMALLEST POSITIVE EPS SUCH THAT 1+EPS <> 1 :   ', EPS:14, ' = ', IBETA, ' ** -', ABS(MACHEP));
  351.    WRITELN ('SMALLEST POSITIVE EPS SUCH THAT 1-EPS <> 1 :   ', EPSNEG:14, ' = ', IBETA, ' ** -', ABS(NEGEP));
  352.    WRITELN ('SMALLEST POS. NORMALIZED FLOATING-POINT NUMBER:', XMIN:14, ' = ', IBETA, ' ** -', ABS(MINEXP)) ;
  353.    WRITELN ('LARGEST FLOATING-POINT NUMBER:                 ', XMAX:14, ' = ', IBETA, ' ** +', ABS(MAXEXP), ' - 1');
  354.    WRITELN;
  355.    WRITE   ('FLOATING-POINT ARITHMETIC ');
  356.    CASE IRND MOD 3 OF
  357.        0: WRITELN ('CHOPS');
  358.        1: WRITELN ('ROUNDS, BUT NOT IEEE STYLE,');
  359.        2: WRITELN ('ROUNDS IEEE STYLE');
  360.    END;
  361.    IF IRND DIV 3 = 1 THEN
  362.       WRITELN ('AND SUPPORTS GRADUAL UNDERFLOW')
  363.    ELSE
  364.       WRITELN ('AND DOES NOT SUPPORT GRADUAL UNDERFLOW');
  365. END; { PrintParam }
  366.  
  367. END. { MachArit }
  368.